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1.  INTRODUCTION 


The  following  sections  present  work  completed  under  the  UNM/AFRL  research  grant  FA9453- 
11-1-0276  since  May  2013.  During  this  time  period,  the  Generalized  Boltzmann  Fokker-Plank 
(GBFP)  advanced  physics  models  were  expanded  to  include  both  discrete  and  hybrid  differential 
cross  section  (DCS)  models.  In  addition,  the  analog  models  that  are  the  basis  of  the  GBFP 
physics  were  expanded  to  include  the  following  DCSs:  the  screened  Rutherford  elastic  scattering 
DCS,  the  partial  wave  expansion  elastic  scattering  DCS,  the  Rutherford  inelastic  DCS,  and  the 
Moller  inelastic  DCS.  Tools  for  tallying  quantities  like  dose,  charge  deposition,  transmission  and 
reflection,  angular  distributions,  energy  spectra,  and  lateral  and  longitudinal  distributions  were 
implemented  and  tested.  Finally,  a  test  suite  including  a  wide  variety  of  problems  was 
completed.  The  remainder  of  the  paper  will  cover  background  material,  discuss  the  methods 
utilized,  and  present  results  obtained  from  the  test  suite. 

2.  BACKGROUND 

Solution  to  the  Boltzmann  transport  equation  for  electrons  requires  information  about  the 
interactions  that  an  electron  can  undergo.  This  information  is  captured  by  the  DCSs  and 
describes  the  probability  that  a  particle  will  scatter  through  some  angle  or  lose  some  energy  to 
the  medium.  In  addition,  the  total  cross  section  is  obtained  from  the  DCS  and  characterizes  the 
length  scale  of  the  physics.  For  electrons,  the  mean  free  path  (MFP)  or  the  inverse  of  the  total 
cross  section  is  extremely  small,  so  particles  suffer  thousands  of  collisions  while  slowing  down. 
Also,  the  DCSs  are  peaked  about  small  changes  in  the  state  of  the  electron.  As  a  result,  analog 
Monte  Carlo  is  extremely  computationally  inefficient  for  electrons  above  a  few  hundred  keV. 

Approximate  methods  were  developed  to  alleviate  the  computational  effort  required  for  analog 
Monte  Carlo  simulation  of  electrons.  However,  methods  like  condensed  history  (CH)  introduce 
error  on  the  order  of  the  fixed  distance  traveled  by  the  electron  between  collisions  or  the  step  and 
CH  requires  special  boundary  crossing  algorithms  to  mitigate  additional  error  incurred  at 
material  interfaces.  The  GBFP  method  is  an  approximate  method,  but  unlike  CH,  the  GBFP 
method  is  a  single-event  Monte  Carlo  method.  That  is,  the  distance  between  collisions  or  the  step 
is  exponentially  distributed,  so  no  boundary  crossing  algorithms  are  required.  In  addition,  the 
GBFP  method  introduces  moment-preserving,  approximate  DCSs  with  systematically 
controllable  accuracy  by  preservation  of  more  moments. 

Like  analog  Monte  Carlo,  the  GBFP  physics  models  require  a  total  cross  section  and  a  DCS  for 
each  interaction  simulated.  The  cross  sections  are  used  to  sample  distance  to  collision,  scattering 
angle,  and  energy-loss.  Not  only  are  the  approximate  DCSs  critical  to  the  GBFP  method,  but  also 
the  analog  DCSs  are  required  to  construct  the  approximate  DCSs  using  moments  of  the  analog 
DCSs.  Therefore,  some  attention  is  given  to  the  DCSs  currently  implemented. 

The  predominant  interactions  that  are  considered  in  this  work  are  elastic  scattering  with  target 
nuclei  and  inelastic  scattering  with  orbital  electrons.  Elastic  scattering  is  a  Coulombic  interaction 
between  primary  electrons  and  target  nuclei  that  can  result  in  a  change  in  the  primary  electron’s 
direction  (energy  loss  due  to  elastic  scattering  is  negligible).  Inelastic  scattering  is  a  Coulombic 
interaction  between  the  primary  electrons  and  the  atomic  electrons  where  energy  from  the 
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primary  electrons  is  transferred  to  the  atomic  electrons  (change  in  direction  due  to  inelastic 
scattering  of  the  primary  electrons  is  neglected).  Currently,  elastic  scattering  models 
implemented  include  the  screened  Rutherford  DCS  [1]  and  the  partial-wave  expansion  DCS  [2], 
The  screened  Rutherford  DCS  is  a  widely  applicable,  but  simplified  model.  The  benefit  of  using 
the  screened  Rutherford  cross  section  is  that  it  is  in  a  continuous  form  or 


C(E) 

P-ft+2  rKE)]2  ’ 


(2.1) 


where  p0  =  cos(0o)  is  the  deflection  cosine,  rj{E)  is  the  screening  parameter,  andC(/t)is  the 
material  constant.  The  screening  parameter  is  given  by 


Z 


2/3 


V(E)  = 


i.i3+3.76(^/37)2(r+i)2rI(r+2)1 


(0.885-137)27Xr+2) 


(2.2) 


where  Z  is  the  atomic  number  of  the  target  nucleus  and  T  =  E  /  m0c2is  the  particle  energy  per 
rest  mass  energy.  For  electrons,  m^c  =  0.510998910MeV.  The  material  constant  is  given  by 


2^0V^z(z+i)  (r+i)2 

y  )  a  t\t+2)2’ 

where  A  is  the  atomic  mass,  r0  =  2.8179403267x10  13  cm  is  the  classical  electron  radius,  p  is 
the  material  density,  and  N. is  Avogadro’s  number.  In  contrast  to  the  screened  Rutherford  DCS, 
the  partial-wave  DCSs  are  given  by 


(2.4) 


where 


=  -t-  52  H/+ 1  )l  exp( 2«TA1 )  - 1  ] + /|  exp<  2/V>;  , )  - 1 1 ; 


(2.5) 


and 


M = ^Xfexp(2/d; , ) - cxp(2/d;  ,)]/)’ (//).  (2.6) 

2  ikI=0 

It  is  beyond  the  scope  of  this  report  to  discuss  the  process  of  calculating  the  partial-wave  DCS 
(see  ELSEPA  [2]);  however,  it  is  important  to  note  that  there  is  no  analytical  form  of  eq.  2.4. 
Therefore,  computer  codes  are  used  to  generate  partial-wave  DCS  libraries.  The  motivation  for 
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using  the  partial-wave  DCS  is  because  it  is  the  most  accurate  representation  of  elastic  scattering 
currently  available.  It  also  provides  an  additional  application  to  demonstrate  the  effectiveness 
and  flexibility  of  the  GBFP  method.  That  is,  the  GBFP  method  does  not  depend  on  the  form  of 
the  DCS,  rather  the  moments.  Therefore,  DCSs  can  be  continuous  or  discrete  so  long  as  moments 
of  the  DCS  can  be  accurately  calculated. 

There  are  two  analog  inelastic  scattering  models  currently  implemented  in  the  Geant4  [4]  test 
code.  The  inelastic  DCS  models  currently  implemented  includes  the  Rutherford  DCS  [3]  given 
by 


Ze(E,Q) 


K(E ) 

e2  ’ 


where  the  material  constant,  K(E),  is 


K(E)  = 


2nr2pNA7, 
Ap 2  ’ 


(2.7) 


(2.8) 


P=  vie,  is  the  particle  velocity  per  c  (c  being  the  speed  of  light),  and  Q  is  the  energy 
transferred  from  the  incident  particle  to  the  atomic  electron.  The  more  accurate  Moller  [3]  DCS 
was  also  implemented  and  is  given  by 


7e(E,Q)=K(E) 


1 


1 


1 


Q 2  {E-QY  (E +m0c2Y 


m0c2(2E  +  m0c2) 


(2.9) 


Given  the  aforementioned  DCS  models,  the  goal  of  the  Monte  Carlo  method  is  to  solve  the 
Boltzmann  transport  equation  for  the  selected  physics  by  simulating  the  particle  behavior 
according  the  DCSs.  For  example,  the  Boltzmann  transport  equation  for  electrons  that  can 
undergo  the  interactions  detailed  above  is 

Q •  V ip(E, Q) + 7t(E) v(E,Q)  =  f  Q)i^(E,Q')+  V dE'Y.  _(£'->£M£',n).  (2.10) 

J  An  J  0  e 

In  the  eq.  2.10  the  spatial  dependence  is  implied  for  compactness,  but  in  the  following  definitions 
the  spatial  dependence  is  explicit: 

yAr,  E,  £2)  Angular  flux 

Hn(r,E,/l0)  Macroscopic  elastic  scattering  cross  section 
UrXQ)  Macroscopic  inelastic  scattering  cross  section 
l^r^E)  Total  macroscopic  interaction  cross  section 

where  r  is  position  vector,  E  is  the  particle  energy,  and  Q  the  direction  of  the  particle.  The 
prime  notation  on  E  and  Q  corresponds  to  the  particles  energy  or  direction  before  a  collision 
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occurs.  The  left-hand  side  of  eq.  2.10  corresponds  to  particles  leaving  the  phase  space  or  particle 
losses  and  the  right-hand  side  corresponds  to  particles  entering  the  phase  space  or  particle  gains. 
The  first  loss  term, 


Q-  V  i fHr,E,Q)  ^Leakage, 

accounts  for  particles  leaking  or  leaving  the  differential  volume.  The  second  loss  term 

Z,(£)y<r,£,Q)  =Outscatter, 

accounts  for  outscatter  or  particles  leaving  the  energy  or  angular  space  by  changing  direction  or 
energy.  The  first  gain  term, 

f  (£,  Q  '•  O)  O  ')=  Elastic  Inscatter, 

J  4^- 

accounts  for  particles  scattering  into  the  angular  phase  space.  The  second  gain  term, 

J  dE"Le  (£’’->■£’) yj{E\Ci)= Inelastic  Inscatter, 

accounts  for  particles  scattering  into  the  energy  phase  space.  According  to  eq.  2.10  and  a 
corresponding  boundary  condition  or  fixed  source,  a  source  particle  is  created  with  a  position, 
energy,  and  direction.  Next,  the  particle  streams  or  moves  a  distance  to  a  collision  site.  This 
process  is  sampled  according  to  the  following  probability  distribution  function 

P{s)ds  -  E/  exp(-X,s)  ds ,  (2.11) 


where  5  is  the  distance  to  collision.  Given  a  distance  to  collision,  the  particle  is  moved  to  a 
collision  site.  At  the  collision  site,  the  collision  type  is  randomly  sampled  according  to  the 
probability  of  the  respective  collision  type.  The  probability  of  an  elastic  scatter  is  the  ratio  of  the 
total  elastic  scatter  cross  section,  X  0,  to  the  total  interaction  cross  section,  X,  =  Xn0  +  Xp0,  or 


P  = 


(2.12) 


and  the  probability  of  an  inelastic  scatter  is  the  ratio  of  the  total  inelastic  scatter  cross  section, 
X^0,  to  the  total  interaction  cross  section,  X/?  or 


P  = 


2, 


(2.13) 
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Given  the  collision  type,  the  collision  outcome  is  determined  by  sampling  the  corresponding 
DCS.  This  is  accomplished  by  direct  or  rejection  sampling  techniques.  For  example,  the  screened 
Rutherford  elastic  scattering  DCS  is  invertible,  so  a  direct  technique  is  used  or 


1-%+T} 


(2.14) 


where  g  is  a  random  number  sampled  uniformly  between  zero  and  one.  However,  the  Moller 
DCS  is  not  invertible  so  a  rejection  technique  must  be  used.  Rejection  techniques  are  well 
documented  for  the  Mdller  DCS  (see  Geant4  Physics  Reference  Manual  [5])  and  not  summarized 
here  because  it  is  an  involved  process.  Regardless  of  the  technique,  the  outcome  is  the  same. 

That  is,  an  energy  loss  is  sampled  according  to  the  physics  governing  the  collision  type  or  the 
DCS.  This  process  of  sampling  distance  to  collision,  collision  type,  and  collision  outcome  is 
repeated  for  each  source  particle  until  all  source  particles  are  complete.  During  the  simulation 
process,  tallies  of  interest  are  accumulated.  As  previously  mentioned,  we  are  interested  in 
quantities  like  dose,  charge  deposition,  transmission  and  reflection,  angular  distributions,  energy 
spectra,  and  lateral  and  longitudinal  distributions.  To  tally  the  various  quantities  of  interest,  a 
running  sum  is  accumulated  as  the  particles  tranverse  the  medium.  For  example,  in  a  one¬ 
dimensional  slab,  dose  in  MeV  cml  per  source  particle  in  the  ith  cell  is 

g 


D  = 


NpAxt  . 


(2.15) 


where  dUj  is  the  dose  deposited  in  the  ith  cell  by  the  jth  source  particle,  Ards  the  width  of  the  ith 
cell,  and  N  is  the  total  number  of  source  particles. 


This  simple  description  of  the  Monte  Carlo  algorithm  is  valid  for  analog  or  single-event  Monte 
Carlo  where  collision  sites  are  exponentially  distributed.  One  significant  difference  between  the 
GBFP  method  and  other  approximations  is  that  the  GBFP  method  is  a  single-event  method. 
Therefore,  pre-existing  single-event  Monte  Carlo  algorithms  can  be  used,  so  no  code 
modifications  are  necessary  to  use  this  approximation.  That  is  because  the  GBFP  method  simply 
relies  on  approximate  DCSs.  There  is  a  very  subtle  difference  between  the  transport  equation 
solved  using  the  GBFP  method  and  the  analog  transport  equation  in  eq.  2.10.  After  substituting 
the  approximate  cross  sections  the  transport  equation  becomes 

Q •  V ip(E, Q) + *KE,Q)  =  f  dQ'tn(E,Q'QME,Q')+  V dE't  _{E' ^ E)yKE' ,&) .  ( 2.16) 

J  J 0  e 

where  2  implies  that  the  cross  section  is  approximate.  Notice  that  the  form  of  operators  in  eq. 
2.10  and  eq.  2.16  are  the  same.  For  this  reason,  many  of  the  issues  encountered  when  using  CH 
methods,  like  boundary  crossings  and  step-length  limitations,  are  avoided  when  using  the  GBFP 
method.  Unlike  CH,  the  transport  equation  remains  unchanged.  The  only  question  that  remains  is 
how  does  one  construct  a  GBFP  DCS  such  that  a  solution  eq.  2.16  is  a  good  approximation  of  the 
analog  solution. 


5 

Approved  for  public  release;  distribution  is  unlimited. 


To  do  so,  Lewis  theory  [6]  is  utilized.  Lewis  theory  provides  a  relationship  between  Legendre 
moments  of  the  elastic  scattering  DCS  and  space-angle  moments  of  the  solution.  That  is,  if  an 
approximate  DCS  preserves  a  finite  number  of  Legendre  moments,  L,  of  the  analog  DCS,  then  it 
is  possible  to  preserve  p+q=L  space-angle  moments  of  the  analog  solution,  where 

Ml  =2x\ P^SE^dp,  (2.17) 

-1 

are  the  Legendre  moments  of  the  analog  elastic  scattering  DCS.  The  space-angle  moments  of  the 
solution  are 

(2-18) 

-CO  -1 

Lewis  theory  is  only  valid  in  an  infinite  medium,  but  it  does  emphasize  the  relationship  between 
moments  of  the  DCS  and  moments  of  the  solution.  Therefore,  the  GBFP  DCSs  are  constructed 
such  that  they  preserve  moments  of  the  analog  DCSs.  To  do  so,  we  write  down  a  system  of 
equations, 


Mj  =MP  1  =  12,...,L  (2.19) 

where 

M,  =2x\ P^lSE^dp.  (2.20) 

-1 

In  eq.  2.19,  L  is  determined  by  the  number  of  unknown  parameters  used  to  construct  the  GBFP 
DCS.  For  both  the  discrete  DCS  and  the  hybrid  DCS,  there  are  a  total  of  27Vpoints  and  weights 
that  must  be  determined.  Thus,  L  —  2AGbr  the  system  in  eq.  2. 19.  As  an  example,  we  take  a 
discrete  DCS  with  the  following  form 


(2.21) 

1,-1 

where  gn  are  the  discrete  points  and  an  are  the  corresponding  weights.  Substituting  eq.  2.21  into 
eq.  2.19  gives 


Ml=a1{E)Pl(g)+...+aN{E)PI(gN),  f  =  l2,-,2 N.  (2.22) 

We  now  have  a  system  of  2 /Vpoints  and  weights  that  we  can  solve  such  that  the  first  moment 
and  higher  order  moments  of  the  discrete  DCS  are  equivalent  to  the  moments  of  the  analog  DCS. 
Therefore,  the  zeroith  moment  or  the  total  cross  section  is  not  preserved.  The  result  is  a  smoother 
cross  section  with  a  longer  mean  free  path  (MFP),  which  are  both  features  of  the  analog  DCS 
that  result  in  substantial  computational  cost  for  analog  Monte  Carlo  simulations.  A  smoother 
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DCS  with  a  longer  MFP  that  preserves  moments  of  the  analog  DCS  is  the  key  to  accurately  and 
efficiently  approximating  eq.  2.10. 

The  goal  of  the  GBFP  method  is  to  utilize  approximate  DCSs  that  are  smoother  and  have  longer 
MFPs.  However,  there  is  a  certain  level  of  detail  that  is  required  to  obtain  these  ideal  DCSs.  The 
following  section  provides  a  thorough  discussion  of  the  process  of  obtaining  various  GBFP 
DCSs. 

3.  METHODS,  PROCEDURES,  AND  ASSUMPTIONS 

The  current  implementation  of  the  GBFP  method  allows  for  two  approximate  DCSs:  discrete  and 
hybrid  DCS  models.  The  discrete  DCS  is  a  set  of  discrete  points  and  weights.  Whereas,  the 
hybrid  DCS  is  a  combination  of  a  smooth  tail  and  a  discrete  peak.  Both  DCSs  require  calculation 
of  discrete  points  and  weights,  so  the  following  sections  discuss  the  general  methods  used  to 
generate  the  DCSs  and  the  details  specific  to  the  discrete  and  hybrid  DCSs.  In  addition,  Geant4 
implementation  details  relevant  to  the  work  completed  to  date  are  presented. 

3.1.  Quadrature  Methods 

In  section  2,  the  system  of  equations  in  eq.  2.21  demonstrates  the  GBFP  concept,  but  that  system 
is  not  directly  invertible  due  to  numerical  instabilities.  Therefore,  we  turn  to  a  method  typically 
used  to  generate  quadrature  sets  as  an  alternative  approach.  The  aim  of  Gaussian  quadrature  is  to 
obtain  the  points,  xN),  and  weights,  (w1,w2,...,wAr),  such  that 

j  *  w(V)JXv)  dfi = X  wJ(xi )  (31) 

F=1 

is  equivalent  for  polynomials,  f(ji),  of  order  2N -  land  less.  To  obtain  the  points  and  weights 
one  must  generate  N  polynomials  that  are  orthogonal  with  respect  to  the  weight  function,  w(p). 
The  polynomials  are  then  used  to  obtain  the  points  and  weights.  That  is,  the  zeros  of  the  order  N 
polynomial  are  the  points  and  the  weights  are  given  by 


_  (/Vil/V-.) 

Pn  i(-0pV(-0  . 


(3.2) 


To  generate  orthogonal  polynomials,  the  recurrence  coefficients,  (tf-,/7),  satisfying 

PJ+ i(x)  =  (x-a^pjix)- f}jApH(x) 


(3.3) 


are  required.  In  the  context  of  this  discussion,  there  are  two  types  of  orthogonal  polynomials, 
classical  and  non-classical.  Classical  polynomials  are  orthogonal  with  respect  to  classical  weight 
functions  (examples  include  Legendre,  Jacobi,  Chebyshev,  Laguerre,  Hermite,  and  so  on).  There 
are  typically  well-known  recurrences  and  recurrence  coefficients  for  classical  orthogonal 


7 

Approved  for  public  release;  distribution  is  unlimited. 


polynomials.  That  is,  computation  of  the  recurrence  coefficients  is  unnecessary,  so  generating 
Gauss-Lengendre  quadrature,  Gauss-Jacobi  quadrature,  and  so  forth  is  just  a  matter  of 
determining  the  points  and  weights  from  the  known  recurrence  coefficients.  However,  it  can  be 
advantageous  to  generate  quadrature  sets  corresponding  to  non-classical  weight  functions. 
Therefore,  an  additional  step  is  required  where  the  recurrence  coefficients  are  obtained  for  the 
polynomials  that  are  orthogonal  with  respect  to  the  non-classical  weight  function. 

It  is  possible  to  obtain  the  recurrence  coefficients  for  non-classical  orthogonal  polynomials  using 
the  modified  Chebyshev  Algorithm  (MCA).  The  MCA  is  derived  by  Gautschi  [7]  and  is  a 
mapping  procedure  for  mapping  the  recurrence  coefficients  of  classical  orthogonal  polynomials 
to  the  recurrence  coefficients  non-classical  orthogonal  polynomials  with  modified  moments 
given  by 


v,  =  J  71  -(x)w(x)dx ,  (3.4) 

J  J  a  J 

where  7Tj(x)  is  a  classical  orthogonal  polynomial  of  order  j .  Given  the  modified  moments  and 

the  recurrence  coefficients  of  classical  orthogonal  polynomials,  the  MCA  returns  the  recurrence 
coefficients  of  non-classical  orthogonal  polynomials  that  are  used  to  generate  the  corresponding 
quadrature  set.  The  MCA  is  typically  used  when  integrating  specialized  functions  that  cannot  be 
factored  into  a  classical  weight  function  and  a  polynomial.  However,  the  system  encountered 
when  generating  a  GBFP  DCS  can  be  recast  into  a  Gaussian  quadrature  problem  for  a  non- 
classical  weight  function. 

Before  moving  to  the  details  of  applying  the  MCA  to  the  GBFP  method,  we  will  write  down  the 
final  step  in  generating  the  points  and  weights.  There  are  a  number  of  methods  for  obtaining  the 
discrete  points  and  weights,  but  the  method  used  in  this  work  is  due  to  Golub  and  Welsch  [8], 
Wilf  recognized  that  the  Gaussian  quadrature  system  could  be  written  in  terms  of  an  eigenvalue 
problem.  The  important  result  is  that  the  zeros  of  the  order  N  polynomial  are  equivalent  to  the 
eigenvalues  of  the  Jacobi  matrix 


a0 

M 

0 

... 

0 

M 

a. 

M 

0 

M 

M 

• 

M 

a3 

0 

v'A 

0 

... 

0  x/a7 

aN 

where  )  are  the  recurrence  coefficients  of  the  orthogonal  polynomials  (classical  or  non- 

classical).  The  first  entry  of  the  Ith  eigenvector,  q.,  squared  corresponds  to  the  z*  weight  or 
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(3.6) 


Finally,  it  is  important  to  note  that  the  GBFP  method  requires  Gauss-Radua  quadrature  where 
one  of  the  points  is  constrained  to  be  the  upper  or  lower  bound  of  the  integral  in  eq.  3.1.  Given 
this  constraint,  the  process  of  obtaining  the  correct  Jacobi  matrix  is  outlined  by  Golub  [9],  The 
important  results  are  summarized  below.  We  wish  to  determine  a  polynomial  such  that 
pN+1(a)=0  or  pN+](b )  =  0,  where  a  and  b  are  the  limits  of  the  integral  in  3.1.  This  implies  that 


(3.7) 


or 


(3.8) 


Now,  the  Jacobi  matrix  is  expanded  from  an  N2  to  an  (/V + 1)2  matrix  with  the  constraint  that 


one  of  the  points,  xP  is  a  or  b.  Given  a  method  for  obtaining  discrete  points  and  weights 

corresponding  to  a  non-classical  weight  function,  it  is  possible  to  recast  the  GBFP  cross  section 
construction  problem  into  a  similar  quadrature  problem  of  a  non-classical  weight  function. 

3.2.  GBFP  Discrete  Cross-section  Generation 

We  can  rewrite  eq.  3.1  in  the  following  form: 


(3.9) 


n= 1 


We  would  like  to  satisfy  this  equation  for  polynomials  of  order  2N— 1  or  less.  Specifically,  we 
would  like  to  satisfy  the  above  equation  for  f(ji)=Pl(jj),  where  Pt(ji)  is  a  Legendre 

polynomial.  Since  eq.  3.9  holds  for  all  polynomials,  the  same  will  be  true  if  f(/i)  =  PI(ji)  or 


(3.10) 


Given  the  discussion  in  the  previous  section,  it  should  be  clear  that  eq.  3.9  and  3.10  can  be 
satisfied  and  the  process  of  mapping  the  known  recurrence  coefficients  to  those  of  non-classical 
recurrence  coefficients  using  the  MCA  is  a  viable  approach.  These  particular  non-classical 
recurrence  coefficients  correspond  to  polynomials  that  are  orthogonal  with  respect  to  the  DCS 
that  is  used  as  a  weight  function.  After  applying  the  Golub  and  Welsch  algorithm  to  these 
recurrence  coefficients,  the  corresponding  points  and  weights  will  satisfy  eq.  3.10  for 
l  =  0,...,2N- 1.  That  is,  the  first  2  N  moments  of  the  DCS  are  preserved.  However,  as  discussed, 
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the  goal  is  to  specifically  avoid  preservation  of  the  zeroith  moment,  because  the  zeroith  moment 
gives  the  total  cross  section  or  the  MFP  causing  computational  inefficiency.  This  is  the 
motivation  for  using  the  Radau  method  discussed  above.  By  requiring  one  discrete  point  is  equal 
to  unity,  we  generate  a  deflection  cosine  that  corresponds  to  a  non-interaction  ( —  I  =  cos(O)). 

A  scattering  angle  of  zero  is  equivalent  to  not  interacting.  Since  the  analog  DCS  is  peaked  about 
non-interaction,  we  seek  to  remove  this  feature  from  the  GBFP  DCSs.  After  generating  the 
discrete  points  and  weights  only  the  points  not  equal  to  unity  are  kept. 

To  make  this  discussion  clear,  we  will  provide  an  example  of  the  process  of  constructing  a  GBFP 
elastic  scatting  DCS.  First,  we  must  select  a  few  parameters.  We  will  demonstrate  the  process  of 
generating  4  discrete  points  and  weights  for  a  1  MeV  electron  in  Gold  with  the  following 
parameters: 


•  Z=79 

•  A=  196.97  g/mol 

•  p= 19.33  g/cc 

•  /7=le-4 

•  N=4 


To  generate  polynomials  that  are  orthogonal  to  the  screened  Rutherford  DCS  we  need  2N+2 
moments  given  by 


M,  =2ir\ 

-1 


(3.11) 


There  are  several  methods  that  can  be  used  to  generate  these  moments  accurately.  One  is  the 
Spencer  recurrence  [10].  This  recurrence  is  only  good  for  the  screened  Rutherford  DCS  over  the 
full  range.  We  implemented  a  more  general  method  such  that  any  elastic  scattering  DCS  for  the 
full  or  partial  range  is  possible.  This  method  simply  uses  adaptive  quadrature  to  evaluate  the 
integral  in  eq.  3.11  and  gives  the  following  for  the  first  10  moments: 

Table  1:  Moments  of  the  Screened  Rutherford  Analog  DCS 


/ 

(l/cm) 

/ 

(l/cm) 

0 

13222.95284269771 

5 

12993.68379391062 

1 

13200.72491556212 

6 

12920.77698622844 

2 

13164.40260046155 

7 

12841.46869195591 

3 

13116.67842732473 

8 

12756.48220266539 

4 

13059.33200551864 

9 

12666.43977610700 

We  must  now  generate  the  first  2N+\  recurrence  coefficients  for  monic  Legendre  polynomials. 
The  coefficients  are  given  by 
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where  /?0  is  not  used  and  can  be  set  to  an  arbitrary  value.  Given  Mv  ,  and  we  can  now 
generate  the  non-classical  recurrence  coefficients  a}  and  bj .  The  result  is 


Table  2:  Non-classical  Recurrence  Coefficients  for  Example  Problem 


j 

aj 

b, 

0 

0.9983189891546911 

arbitrary 

1 

-0.004950428124079198 

0.0004072459551613386 

2 

-0.002528000016159826 

0.3308498723933296 

3 

-0.001661056221544444 

0.2655511190282964 

4 

-0.001207307704623406 

0.2564056575870433 

The  coefficients  in  table  2  are  passed  to  the  Radau  algorithm  and  the  discrete  points  and  weights 
are  generated.  In  table  3  and  table  4,  the  discrete  points  and  weights  are  given  before  and  after 
regularizing  or  canceling  the  non-interaction  point  respectively. 

Table  3:  Example  Discrete  Cross-Section  Before  Regularizing 


i 

C 

l 

-0.8253336745219344 

0.3542588039339951 

2 

-0.1926588357286153 

1.46098457679131 

3 

0.5530858576077967 

8.8112540880515 

4 

0.974594146260808 

625.8776637056218 

5 

0.9999999999999994 

12586.44868152331 

Table  4:  Example  Discrete  Cross-Section  After  Regularizing 


/ 

C 

<E)n 

1 

-0.8253336745219344 

0.3542588039339951 

2 

-0.1926588357286153 

1.46098457679131 

3 

0.5530858576077967 

8.8112540880515 

4 

0.974594146260808 

625.8776637056218 
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The  definition  of  the  total  cross  section  is  the  sum  of  the  weights  or  13222.9528  (1/cm)  prior  to 
regularizing.  After  regularizing,  the  total  cross  section  is  636.5042  (1/cm).  This  means  that  the 
MFP  associated  with  the  regularized  cross  section  is  about  21  times  longer  than  the  analog  MFP 
or  on  average  a  1  MeV  analog  electron  in  gold  suffers  2 1  times  more  collisions  than  a  GBFP 
electron. 

So  far,  we  have  discussed  the  elastic  scattering  DCS  construction  using  quadrature  methods. 
However,  these  methods  are  only  valid  over  [-1,1].  Therefore,  construction  of  a  GBFP  inelastic 
scattering  DCS  requires  additional  steps.  It  is  possible  to  determine  a  linear  mapping  between  Q, 
the  dependent  variable  for  inelastic  scattering,  and  ju,  the  dependent  variable  for  elastic 
scattering.  That  is,  the  following  is  true  for  the  elastic  and  inelastic  scattering  DCSs: 


UmdQ= 27rln(E,MW- 


(3.13) 


Given  a  mapping,  one  can  compute  Legendre  moments  of  the  mapped  inelastic  scattering  DCS. 
These  moments  are  then  used  to  generate  points  and  weights  on  [-1,1]  and  then  the  points  and 
weights  are  mapped  back  to  Q  space.  This  process  is  outlined  below. 

First,  we  must  determine  a  mapping.  If  we  assume 


(3.14) 


and  require  that  Q(ji — 1)  =  Qmax  and  Q(ji  =  - 1)  =  0,  we  determine  m  and  b  such  that 


(3.15) 


Given  this  mapping,  we  can  now  obtain  a  relationship  between  energy- loss  moments,  Qn,  and 
Legendre  moments  of  a  pseudo  elastic  scattering  DCS  or  Mv  We  would  like  to  rewrite  the  right 
hand  side  of  the  following  equation  in  terms  of  energy  loss  moments 


M,  =  2*\\  P,(ji)I.n(E,fi)dfi. 


(3.16) 


From  our  mapping,  we  know  that  (1—  fi)=  Vo  Q  ancl we  can  rewrite  P{(ji)  as 


(3.17) 


where 


c 


n 


-LfW-D-V-l)] 

2  "!  i=o 


(3.18) 
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Equation  3.16  becomes 


M l  =  t-yj-  J  ^(1-z^y  £„(£, ii)dfi . 

j= i  J' 

Then  we  substitute  the  expression  for  (1— pi) and  Hn(E,pi)dpi  and  change  the  limits  of 
integration,  which  gives 


(3.19) 


M,  =  \<’~Q,'£,(k,Q)dQ 

j=l  J' 


(3.20) 


or 


^iyS^-^ojQr  <3-2» 

j= i  J- 

where 

=  (3.22) 

are  the  energy-loss  moments  of  the  inelastic  DCS. 

The  procedure  for  generating  a  GBFP  inelastic  DCS  is  the  same  as  for  the  elastic  DCS  except  we 
start  with  energy-loss  moments.  After  mapping  the  moments,  we  can  use  the  same  procedure 
outlined  earlier  in  this  section.  This  will  result  in  points  that  are  on  [-1,1]  and  we  simply  map 
those  points  back  to  Q  space  using  eq.  3.15  to  complete  the  process. 

3.3.  GBFP  HYBRID  CROSS-SECTION  GENERATION 

Generation  of  the  hybrid  DCS  requires  the  same  procedure  discussed  above,  but  the  interval  over 
which  the  discrete  points  are  generated  is  restricted  to  1/^1]  for  elastic  scattering  and 

[£L>£U  for  inelastic  scattering.  The  remaining  interval,  [-1,/i  J  and  [Qcul,QmJ,  is 

represented  exactly  by  the  respective  analog  DCS.  The  first  parameter  that  must  be  determined  is 
the  cut  value.  The  cut  value  is  driven  by  the  underlying  principle  of  the  hybrid  DCS.  That  is,  by 
representing  the  catastrophic  collisions  (those  resulting  in  large-angle  deflection  and  large  energy 
loss)  with  the  analog  DCS,  discrete  artifacts  are  mitigated.  Therefore,  selection  of  the  cut  should 
relate  to  how  much  of  the  GBFP  DCS  is  represented  by  the  analog  DCS.  To  do  so,  we  select  the 
cut  such  that  total  cross  section  corresponding  to  the  smooth  portion  of  the  GBFP  DCS,  Sf ,  is 
some  fraction  of  the  analog  total  cross  section  or 
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X;VQ=^ 

'  N  . 


(3.23) 


We  then  solve  for  the  cut  by  evaluating  or,  for  screened  Rutherford  elastic  scattering, 


fft*  C(£)  X,(£) 

J  1  [1-//+2t7(£)]2  “ 


(3.24) 


This  gives 


^=1+2^- 


N2t](l+ij) 

l+Nr/ 


(3.25) 


If  Nis  unity,  then  the  hybrid  DCS  goes  over  to  the  analog  DCS.  If  N  is  large,  fl^  — 1  and  the 

hybrid  DCS  goes  over  to  discrete  DCS.  The  cut  value  can  be  used  to  control  the  speed  and 
accuracy  of  the  hybrid  model  and  must  be  determined  by  the  applicable  analog  DCS. 


Given  a  cut  value,  the  moments  used  to  obtain  discrete  points  and  weights  for  elastic  scattering 
become 


V,=  f  (3.26) 

Mad 

and  for  inelastic  scattering 

Q,=  \<^Q”Z.(.E,Q)dQ.  (3.27) 


The  procedure  for  generating  elastic  scattering  DCSs  must  change  to  incorporate  a  mapping  from 
1/^1]  to  [-U]  before  utilizing  the  MCA.  A  process  similar  to  what  was  required  for  inelastic 

scattering  is  applied.  Rather  than  going  through  the  derivation  of  the  mapping,  we  will  simply 
provide  the  result.  First,  the  mapping  for  the  moments  is 


vi=1 X 


(-i  r 


n= 0 


n\ 


-ft.  /  Jfc=0 


(3.28) 


where  Mk  is  given  in  eq.  3.26,  cn  is  given  in  eq.  3.18,  and  b"k  is 


\«+m+ 1 


m  0 


m\  n+m+ 1 


If  fi  g  [/G^,l]  and  ji  e  [—1,1],  the  linear  mapping  from  fi  to  fi  is 


(3.29) 
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(3.30) 


The  procedure  for  generating  discrete  inelastic  scattering  DCSs  remains  the  same  except  that  the 
mapping  changes  to 


0C«)=^<  1-/0 

and  the  moments  supplied  to  the  MCA  are  mapped  from  the  moments  in  eq.  3.27. 


(3.31) 
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3.4.  IMPLEMENTATION  OF  THE  GBFP  MODELS  IN  GEANT4  TEST 
CODE 

Geant4  is  a  toolkit  for  the  simulation  of  the  passage  of  particles  through  matter.  Its  areas  of 
application  include  high  energy  nuclear  and  accelerator  physics,  medical  physics,  and  space 
physics  [4],  In  general,  use  of  Geant4  entails  developing  an  application  with  at  least  three 
mandatory  user-defined  classes  that  are  inherited  from  the  base  classes  named 
G4VDetectorConstruction,  G4VPhysicsList,  and,  G4VPrimaryGenerator Action.  Creation  of  a 
G4VDetectorConstruction  inheriting  classes  is  used  to  set  the  problem  geometry.  In  this 
inheriting  class,  the  dimensions  and  material  properties  for  the  model  are  set.  A  G4VPhysicsList 
inheriting  class  is  used  to  set  the  various  physics  processes  that  are  applicable  for  the  particles 
required  by  the  simulation.  In  this  inheriting  class,  the  processes  and  the  associated  models  are 
created  and  applied  to  the  appropriate  particles.  Finally,  creation  of  a 

G4VPrimaryGeneratorAction  inheriting  class  is  used  to  define  the  particle  source.  The  position, 
energy,  particle  type,  and  direction  are  set  in  this  class  (distributions  for  each  of  these  parameters 
can  be  used  as  well).  It  is  at  the  level  of  processes  and  models  where  the  details  of  the 
implementation  of  the  GBFP  method  are  of  interest.  A  process  and  model  was  written  for  each 
GBFP  DCS  (Discrete  and  Hybrid)  and  for  each  interaction  (elastic  and  inelastic).  In  addition,  a 
GBFP  DCS  generation  tool  was  written  and  used  by  the  various  models. 

The  physics  processes  determine  when  and  how  to  apply  the  physics  (at  the  beginning  of  the 
step,  along  the  step,  or  at  the  end  of  a  step).  For  example,  all  of  the  physics  processes  created  for 
the  test  code  are  single  event,  so  the  processes  are  applied  at  the  end  of  a  step.  Depending  on  the 
type  of  process  (elastic  or  inelastic)  the  process  samples  an  angular  deflection  or  energy  loss  and 
then  updates  the  state  of  the  particle  accordingly.  The  physics  models  are  used  to  define  how  to 
sample  the  collision  outcome  and  the  total  cross  section  for  that  particular  interaction.  In 
addition,  the  data  or  the  GBFP  DCSs  are  prepared  in  the  physics  models  before  the  run  begins.  A 
GBFP  DCS  and  total  cross  section  is  generated  for  each  element  in  the  model  at  the  time  the 
physics  models  are  initialized.  However,  one  could  generate  a  GBFP  DCS  database  rather  than 
preparing  the  data  at  the  time  of  the  run. 

The  discrete  GBFP  models  are  contained  in  a  single  class.  However,  the  hybrid  GBFP  model  is  a 
combination  of  a  soft  collision  model  and  a  hard  collision  model.  The  soft  collision  model 
corresponds  to  the  discrete  portion  of  the  DCS  and  has  the  same  structure  as  the  discrete  GBFP 
models.  The  hard  collision  model  corresponds  to  the  smooth  portion  of  the  DCS  and  has  the 
same  structure  as  the  analog  DCS  models. 

Each  of  the  following  GBFP  physics  model  and  process  classes  are  currently  implemented: 

•  AnalogScreenedRutherfordElasticScatteringModel 

•  AnalogScreenedRutherfordElasticScatteringProcess 

•  AnalogMottElasticScatteringModel 

•  AnalogMottElasticScatteringProcess 

•  AnalogRutherfordlnelasticScatteringModel 

•  AnalogRutherfordlnelasticScatteringProcess 

•  AnalogMollerlnelasticScatteringModel 
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•  AnalogMollerlnelasticScatteringProcess 

•  DiscreteScreenedRutherfordElasticScatteringModel 

•  DiscreteScreenedRutherfordElasticScatteringProcess 

•  DiscreteMottElasticScatteringModel 

•  DiscreteMottElasticScatteringProcess 

•  DiscreteRutherfordlnelasticScatteringModel 

•  DiscreteRutherfordlnelasticScatteringProcess 

•  DiscreteMollerlnelasticScatteringModel 

•  DiscreteMollerlnelasticScatteringProcess 

•  HybridSoftScreenedRutherfordElasticScatteringModel 

•  HybridSoftScreenedRutherfordElasticScatteringProcess 

•  HybridSoftMottElasticScatteringModel 

•  HybridSoftMottElasticScatteringProcess 

•  HybridSoftRutherfordlnelasticScatteringModel 

•  HybridSoftRutherfordlnelasticScatteringProcess 

•  HybridSoftMollerlnelasticScatteringModel 

•  HybridSoftMollerlnelasticScatteringProcess 

•  HybridHardScreenedRutherfordElasticScatteringModel 

•  HybridHardScreenedRutherfordElasticScatteringProcess 

•  HybridHardMottElasticScatteringModel 

•  HybridHardMottElasticScatteringProcess 

•  HybridHardRutherfordlnelasticScatteringModel 

•  HybridHardRutherfordlnelasticScatteringProcess 

•  HybridHardMdllerlnelasticScatteringModel 

•  HybridHardMdllerlnelasticScatteringProcess 
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4.  RESULTS  AND  DISCUSSION 


The  following  sections  present  results  obtained  since  the  last  reporting  period.  Namely,  the  more 
interesting  results  from  the  test  suite  are  highlighted.  Each  of  the  following  items  are  discussed  in 
remaining  subsections 

•  CEASE-like  problem 

•  Angular  deflection  distributions, 

•  Energy-loss  spectra, 

•  Dose-depth  curves, 

•  Two-dimension  dose  deposition, 

•  Reflection  and  transmission  fractions. 

The  first  item,  “CEASE-like  problem”,  is  a  demonstration  of  the  GBFP  method  in  one- 
dimension  slabs  with  thicknesses  on  the  order  of  the  CEASE  (compact  environmental  anomaly 
sensor)  silicon  detectors.  In  addition,  a  source  with  a  power  law  energy  distribution  was  tested  on 
these  problems,  since  power  law  spectra  are  encountered  in  the  radiation  belts. 

The  remaining  items  are  results  from  an  extensive  test  suite.  The  test  suite  was  used  to  find 
programming  errors  and  to  characterize  the  performance  of  the  GBFP  method  for  various 
problem  parameters  (e.g.  particle  energy,  material  type,  result  type,  and  so  on). 

All  of  the  problems  in  the  test  suite  have  an  associated  length.  For  finite  medium  problems,  the 
associated  length  gives  a  measure  of  the  thickness  of  the  medium.  The  length  tested  is  related  to 
a  condensed  history  step-size.  The  definition  of  the  step-size  used  in  this  paper  and  in  many 
condensed  history  codes  is 


*,  =  !“'  (-%r  dE, 

'Vi 


(4.1) 


where  £\  is  the  initial  energy  of  the  particle  at  the  beginning  of  the  step,  (~d%c)  is  the  stopping 
power,  and 


(4.2) 


is  the  energy  of  the  particle  at  the  end  of  the  step.  Therefore,  a  step-length  is  the  distance  traveled 
while  slowing  down  from  Et  to  £),.  This  metric  is  used  for  presenting  results  because  it 

provides  a  length  scale  that  is  consistent  with  condensed  history.  That  is,  the  GBFP  method  must 
be  effective  at  these  lengths  to  be  considered  a  viable  alternative.  Results  were  generated  for 
slabs  with  thicknesses  of  one  sub-step,  one  step,  and  10  steps  according  to  the  schematic  in  fig. 

1 .  However,  the  quantities  of  interest  vary  from  angular  distributions  to  energy  spectra  and  dose 
deposition  profiles. 
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Figure  1 :  One-Dimension  Slab  Problem  Schematic. 

Since  many  of  the  following  results  were  generated  simultaneously,  we  report  a  single  efficiency 
measure  that  provides  a  sense  of  the  speed-ups  achieved  with  the  GBFP  method.  The  following 
table  gives  the  ratio  of  the  wall  time  required  for  a  GBFP  simulation  to  the  analog  simulation. 
The  tables  indicate  that  there  is  a  dependence  on  both  the  energy  of  the  source  particle  and  the 
level  of  accuracy  of  the  approximation.  For  higher  energy  particles,  the  analog  DCS  is  more 
peaked  and  the  regularization  process  is  much  more  significant.  The  resulting  DCSs  have  longer 
MFPs  for  higher  energy  particles  than  for  lower  energy  particles  relative  to  the  analog  MFP.  The 
other  contributing  factor  impacting  efficiency  is  the  level  of  accuracy.  The  less  accurate 
approximations  are  one  to  two  orders  of  magnitude  more  efficient  that  the  more  accurate 
approximations.  Ultimately,  there  is  a  clear  trade-off  between  efficiency  and  accuracy.  This 
poses  an  interesting  optimization  problem  that  will  be  addressed  in  future  work.  Regardless,  it  is 
clear  that  the  GBFP  method  results  in  efficiency  gains  of  one  to  two  orders  of  magnitude  over 
analog  Monte  Carlo. 

Table  5:  Ratio  of  Wall  Time  for  Analog  Simulation  to  Wall  Time  for  GBFP  Simulation. 


Energy 

(keV) 

Approximation 

1  Angle, 

2  Angles, 

4  Angles, 

8  Angles, 

1  Energy 

1  Energy 

1  Energy 

1  Energy 

10 

0 

0 

0 

0 

50 

2 

2 

1 

1 

100 

4 

3 

2 

2 

500 

22 

12 

6 

4 

1000 

40 

21 

10 

6 

5000 

178 

89 

42 

19 

10000 

381 

196 

90 

41 
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Table  6:  Ratio  of  Wall  Time  for  Analog  Simulation  to  Wall  Time  for  GBFP  Simulation 

Continued. 


Energy  (keV) 

Approximation 

Hybrid  1  Angle 

Vs-  V 
~  Ao 

1  Energy 

Hybrid  1  Angle 

ys  -V 
~  /100 

1  Energy 

Hybrid  1  Angle 

Vs  -  V 

^7  —  /l  (XX) 

1  Energy 

1  Angle, 
Analog 
Energy 

4  Angles 
Analog 
Energy 

10 

0 

0 

0 

0 

0 

50 

1 

2 

2 

2 

1 

100 

2 

3 

4 

5 

2 

500 

2 

5 

13 

17 

6 

1000 

2 

6 

17 

26 

9 

5000 

2 

6 

24 

51 

26 

10000 

2 

6 

25 

60 

38 
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4.1.  CEASE-LIKE  PROBLEM 


Many  of  the  results  obtained  from  the  test  suite  are  for  the  most  challenging  problem  in  electron 
transport:  a  mono-energetic  pencil  beam  incident  on  a  thin  slab.  Pencil  beams  are  difficult  to 
resolve  because  they  are  localized  in  space,  angle,  and  energy.  Whereas,  many  sources 
encountered  in  real  world  problems  are  distributed  in  space,  angle,  and  energy.  The  GBFP 
method  is  actually  more  accurate  when  simulating  distributed  sources.  In  addition,  the  GBFP 
method  is  more  effective  on  thicker  problems  where  the  source  particles  suffer  enough  collisions 
to  spread  out  in  space,  angle,  angle  and  energy.  We  tested  the  GBFP  method  on  slabs  with 
thicknesses  representative  of  the  CEASE  telescope  silicon  detectors  (that  is,  150  microns  for  dft 
and  500  microns  for  dbt  seen  in  fig.  2  [1 1]).  The  CEASE  is  a  small,  lightweight,  and  low  power, 
spacecraft  instrument  designed  to  measure  the  local  space  radiation  environment  and  generate 
warnings  of  the  space  environment  hazards  of  radiation  damage,  dielectric  charging  and  single 
events  effects  [12].  The  150  micron  slab  is  roughly  15  steps  for  a  100-keV  electron. 


Figure  2:  CEASE  Particle  Telescope  Schematic. 


We  look  at  both  energy  spectrum  and  dose  profiles  since  both  of  these  quantities  are  relevant  to 
response  function  generation.  Response  functions  are  necessary  for  interpreting  raw  satellite 
data,  so  it  is  important  that  the  response  functions  are  accurate.  One  method  for  determining 
response  functions  is  simulating  the  response  of  a  detector  using  a  Monte  Carlo  particle  transport 
code  like  Geant4  or  MCNP.  Assuming  an  extremely  accurate  geometric  model  of  the  detector  is 
used,  inaccuracies  can  only  be  a  function  of  the  physics  employed.  For  this  reason,  it  is  important 
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for  any  approximation  to  properly  capture  quantities  like  energy  spectrum  and  dose  while 
remaining  efficient.  We  will  show  that  under  these  conditions,  the  GBFP  method  performs  very 
well. 

The  following  results  were  generated  from  simulations  with  a  mono-directional,  energy- 
dependent  source.  The  source  energy  is  sampled  from  a  power  law  distribution,  since  a  power 
law  is  representative  of  energy  spectra  encountered  in  the  radiation  belts.  The  source  energies 
range  from  1-keV  to  10000-keV.  The  first  result  presented  is  transmitted  energy  spectrum  in  the 
150-micron  slab.  As  seen  in  the  relative  difference  plot,  each  of  the  approximations  are  well 
within  5%  of  the  benchmark. 


E  (MeV) 
(a)  Spectra 


E  (MeV) 

(b)  Relative  Difference 


Figure  3:  Reflected  Energy  Spectrum  (a)  and  Relative  Difference  (b)  for  Electron  Pencil 
Beam  With  Power  Law  Energy  Distribution  in  150  Microns  of  Silicon. 

Next  we  present  one-dimensional  dose  distributions.  As  seen  in  figs.  4a  and  4b,  and  in  5a  and  5b, 
each  of  the  approximations  are  within  a  few  percent  of  the  benchmark. 


Depth  (cm) 


0  0.005  0.01  0.015 

Depth  (cm) 


(a)  Dose  (b)  Relative  Difference 

Figure  4:  Dose  (a)  and  Relative  Difference  (b)  for  Electron  Pencil  Beam  With  Power  Law 

Energy  Distribution  in  150  Microns  of  Silicon. 
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Depth  (cm) 


(a)  Dose  (b)  Relative  Difference 

Figure  5:  Dose  (a)  and  Relative  Difference  (b)  for  Electron  Pencil  Beam  With  Power  Law 

Energy  Distribution  in  500  Microns  of  Silicon. 


The  single-angle,  single-energy  approximation  is  very  efficient  at  nearly  400  times  more 
efficient  than  analog.  Moreover,  this  approximation  is  within  sufficient  levels  of  accuracy 
making  it  a  viable  option  for  dose  calculations  or  response  function  generation. 
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4.2.  ANGULAR  DEFLECTION  DISTRIBUTIONS 


In  this  section,  we  present  the  angular  distributions  of  electrons  that  are  either  transmitted  or 
reflected.  Though  test  cases  were  completed  for  slabs  with  varying  thicknesses,  we  focus  our 
attention  on  the  results  for  slabs  with  thicknesses  of  one  step  and  discuss  the  impact  of  the 
thinner  and  thicker  slabs.  For  these  particular  problem  types,  a  mono-energetic  electron  source 
normally  incident  on  the  left  face  of  a  slab  at  jc  =  Ois  simulated  for  source  energies  ranging  from 
100-keV  to  10000-keV  in  low-Z,  medium-Z,  and  high-Z  materials  (Silicon,  Copper,  Gold). 
However,  we  only  present  the  GBFP  method  in  low-Z  materials  because  low-Z  targets  are  the 
most  challenging.  Both  elastic  scattering  and  inelastic  scattering  are  simulated  using  the  screened 
Rutherford  DCS  and  the  Moller  DCS  respectively,  along  with  GBFP  approximations  of  these 
DCSs.  The  quantities  of  interest  are  the  angular  distributions  of  reflected  and  transmitted 
particles. 

The  benchmark  solutions  were  obtained  using  the  analog  DCSs  and  are  presented  first  to  provide 
a  sense  of  the  physics.  We  begin  with  the  transmitted  angular  distributions.  In  fig. 

6,  the  transmitted  angular  distributions  in  silicon,  copper,  and  gold  for  100-keV  and  10000-keV 
electrons  and  for  a  slab  thickness  of  one  step  are  presented.  The  two  energies  presented  are  the 
upper  and  lower  energies  that  were  tested  and  capture  the  impact  of  the  source  energy  on  the 
angular  distribution.  For  lower  energy  particles,  the  screened  Rutherford  elastic  scattering  DCS 
is  less  peaked  about  forward  scattering.  In  fact,  for  low  enough  energies  the  screened  Rutherford 
elastic  scattering  DCS  is  nearly  isotropic.  Therefore,  the  resulting  transmitted  angular 
distributions  are  effectively  isotropic.  This  is  because  isotropic  scattering  rapidly  spreads  out  the 
mono-directional  beam  in  angle.  For  higher  energy  particles,  the  screened  Rutherford  elastic 
scattering  DCS  is  highly  peaked  and  therefore  spreading  of  the  beam  occurs  slower.  This  is 
shown  clearly  in  fig.  6  as  the  most  peaked  distributions  correspond  to  the  10000-keV  source  and 
the  more  isotropic  distributions  correspond  to  the  100-keV  source. 

The  anisotropy  of  the  angular  distributions  depends  on  the  peakedness  of  the  DCS,  but  also  on 
the  number  of  collisions  suffered  by  the  particles.  In  thin  slabs,  most  particles  escape  before 
suffering  very  many  collisions  and  the  resulting  angular  distributions  are  more  peaked.  In  thick 
slabs  particles  suffer  many  more  collisions  and  the  initial  beam  is  spread  out  in  angle.  However, 
the  size  of  the  slab  is  not  the  only  mechanism  that  leads  to  a  smoother  angular  distribution.  For 
high-Z  materials,  the  effect  nuclear  screening  by  the  atomic  electrons  is  stronger  which  leads  to  a 
larger  average  deflection  cosine.  That  is,  in  high-Z  materials,  electrons  scatter  through  larger 
angles  on  average  than  in  low-Z  materials  causing  additional  spreading  of  the  beam. 
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Figure  6:  Transmitted  Angular  Distributions  for  100-keV  and  10000-keV  Electrons  in 
Silicon,  Copper,  and  Gold  Slabs  with  Thickness  of  One  Step. 

These  benchmarks  indicate  that  the  peakedness  of  an  angular  distribution  is  strongly  dependent 
on  the  source  particles  initial  energy  and  also  dependent  on  the  atomic  number  of  the  target.  It 
can  be  difficult  to  resolve  the  more  peaked  distributions  with  a  discrete  GBFP  DCS.  While  this  is 
true  for  GBFP  DCSs  with  only  a  few  points  and  weights,  it  is  possible  to  resolve  the  transmitted 
and  reflected  angular  distributions  with  a  higher  order  GBFP  DCS.  In  figs.  7a,  7b,  7c,  and  7d 
transmitted  angular  distributions  and  relative  differences  for  100-keV  electrons  in  silicon 
generated  using  the  GBFP  discrete  and  hybrid  DCSs  are  compared  with  the  benchmark.  It  is 
important  to  note  that  angular  distributions  are  differential  quantities,  so  they  are  difficult  to 
resolve  using  any  approximation;  especially,  for  thin  slab  problems  on  the  order  of  a  stepsize  in  a 
low-Z  material.  That  is,  this  is  a  very  strict  test  of  the  GBFP  method  and  for  thicker  slabs  or 
higher-Z  materials  the  GBFP  method  improves  in  accuracy. 

The  single-angle  DCS  results  in  severe  discrete  artifacts.  Though  a  four-angle  DCS  dramatically 
improves  the  result,  one  must  use  eight  points  to  completely  smooth  out  the  artifacts.  The  eight- 
angle  DCS  is  within  5%  of  the  analog  solution.  An  alternative  to  adding  more  angles  to  remove 
artifacts  is  the  hybrid  DCS.  The  hybrid  DCSs  presented  are  represented  by  a  single  point  for  the 
peak  and  the  remainder  of  the  DCS  represented  by  the  smooth  analog  DCS.  We  only  look  at  a 

few  hybrid  DCSs  to  show  the  effect  of  the  cut  value.  The  most  efficient  hybrid  DCS,  Z*  =  X|-1000 , 

results  in  discrete  artifacts.  This  is  because  the  cut  value  is  too  far  from  unity  and  the  discrete 
portion  of  the  DCS  dominates  elastic  scattering.  By  increasing  the  smooth  total  cross  section  by  a 

factor  of  10,  or  Z*  =  ,  the  discrete  artifacts  are  completely  removed  with  exception  of 

scattering  angles  that  are  roughly  zero.  We  can  refine  the  hybrid  DCS  further,  where  Z*  =  %>, 
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but  there  is  only  a  subtle  gain  in  accuracy  for  additional  cost  in  efficiency.  The  impact  of  energy 
dependence  is  captured  in  the  relative  difference  plot.  For  low  energy  particles,  the  addition  of 
more  discrete  energies  is  almost  negligible  for  the  more  accurate  angular  approximations.  There 
is  some  smoothing  of  the  discrete  artifacts  for  the  less  accurate  discrete  angular  approximations, 
but  the  artifacts  remain.  This  indicates  that  a  single-energy  DCS  may  be  sufficient  when 
resolving  angular  distributions  because  the  accuracy  of  the  angular  approximation  is  more 
important.  The  single-energy  DCS  preserves  the  first  two  moments,  or  stopping  power  and 
energy  straggling,  which  seem  to  capture  the  majority  of  the  energy-loss  physics  necessary  for 
these  cases. 


(a)  Comparison  with  Discrete 


(b)  Discrete  Relative  Differences 
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Figure  7:  Comparison  of  GBFP  and  Benchmark  Transmitted  Angular  Distributions  using 
(a)  Discrete  and  (c)  Hybrid  and  the  Relative  Differences  for  (b)  Discrete  and  (d)  Hybrid  for 

100-keV  Electrons  in  Silicon. 


In  figs.  8a,  8b,  8c,  and  8d  transmitted  angular  distributions  for  10000-keV  electrons  in  silicon 
generated  using  the  GBFP  discrete  and  hybrid  DCS  are  compared  with  the  benchmark.  For 
higher  energies,  the  impact  of  discrete  artifacts  is  greater.  Both  a  single-point  and  four-point 
DCS  result  in  artifacts  that  overwhelm  the  angular  distribution.  At  least  an  eight-point  DCS  is 
required  to  reduce  the  effects  of  the  discrete  artifacts.  Some  discrete  artifacts  are  clear  in  the 
green  curve  for  small  scattering  angles.  However,  even  under  this  extreme  test  case,  the  addition 
of  discrete  energies  smooth  out  these  features.  For  the  hybrid  DCS  the  discrete  artifacts  are 
completely  removed  even  with  the  most  relaxed  hybrid  approximation. 
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(a)  Comparison  with  Discrete 


(b)  Discrete  Relative  Differences 
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Figure  8:  Comparison  of  GBFP  and  Benchmark  Transmitted  Angular  Distributions  using 
(a)  Discrete  and  (c)  Hybrid  and  the  Relative  Differences  for  (b)  Discrete  and  (d)  Hybrid  for 

10000-keV  Electrons  in  Silicon. 


Although  the  discrete  DCS  can  result  in  discrete  artifacts,  these  artifacts  are  not  nearly  as 
significant  when  calculating  integral  quantities  like  dose,  charge  deposition,  or  transmission  and 
reflection  (see  following  sections).  If  non-integral  quantities  like  angular  distributions  are 
necessary,  the  hybrid  DCS  is  effective  at  mitigating  the  discrete  artifacts  while  still  achieving 
some  efficiency  gains  over  analog. 
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In  fig.  9,  the  reflected  angular  distributions  in  silicon,  copper,  and  gold  for  various  electron 
energies  and  for  a  slab  thickness  of  one  step  are  presented.  This  figure  shows  a  result  opposite  of 
the  transmitted  angular  distributions.  That  is,  lower  energy  particles  are  more  likely  to  be 
reflected  so  the  magnitude  of  the  lower  energy  curves  is  greater.  However,  there  is  a  consistency 
between  the  transmitted  and  reflected  angular  distributions.  In  both  results,  the  higher  energy 
angular  distributions  are  anisotropic  (tend  to  be  peak  at  10-20  when  transmitted  degrees  or  1 10- 
120  when  reflected).  However,  the  lower  energy  curves  are  nearly  isotropic. 
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Figure  9:  Reflected  Angular  Distributions  for  100-keV  and  10000-keV  Electrons  in  Silicon, 
Copper,  and  Gold  Slabs  with  a  Thickness  of  One  Step. 


We  now  present  results  for  the  reflected  angular  distributions.  In  figs.  10a,  10b,  10c,  and  lOd 
reflected  angular  distributions  for  100-keV  electrons  in  silicon  generated  using  the  GBFP 
discrete  and  hybrid  DCS  are  compared  with  the  benchmark.  Here,  the  discrete  DCS  results  in 

artifacts  for  each  of  the  DCSs  tested.  For  100-keV  sources,  a  hybrid  DCS  with  Yft  —  ^/ooo  is 

sufficient.  It  is  likely  that  reflected  particles  typically  result  from  hard  collisions  that  are  captured 
by  the  tail  of  the  DCS  or  the  smooth  component.  Therefore,  any  soft  collisions  that  are  captured 
by  the  single  discrete  angle  do  not  contribute  significantly  to  the  solution  and  discrete  artifacts 
are  not  present. 
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(a)  Comparison  with  Discrete 


(b)  Discrete  Relative  Differences 


Figure  10:  Comparison  of  GBFP  and  Benchmark  Reflected  Angular  Distributions  using  (a) 
Discrete  and  (c)  Hybrid  and  the  Relative  Differences  for  (b)  Discrete  and  (d)  Hybrid  for 

100-keV  Electrons  in  Silicon. 


In  figs.  11a,  1  lb,  11c,  and  1  Id,  reflected  angular  distributions  for  10000-keV  electrons  in  silicon 
generated  using  GBFP  discrete  and  hybrid  DCSs  are  compared  with  the  benchmark.  Again,  all  of 
the  discrete  DCSs  tested  resulted  in  artifacts.  For  10000-keV  sources,  a  hybrid  DCS  with 
2?  =  ^Xooo  a^so  sufficient.  The  higher  energy  results  are  statistically  noisier  than  the  low  energy 

results  because  very  few  high-energy  particles  are  reflected.  However,  these  results  still  gives  a 
sense  of  the  resulting  reflected  angular  distributions  when  using  a  hybrid  DCS. 
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(a)  Comparison  with  Discrete 


(b)  Discrete  Relative  Differences 


0  (degrees)  0  (degrees) 

(c)  Comparison  with  Hybrid  (d)  Hybrid  Relative  Differences 


Figure  11:  Comparison  of  GBFP  and  Benchmark  Reflected  Angular  Distributions  using  (a) 
Discrete  and  (c)  Hybrid  and  the  Relative  Differences  for  (b)  Discrete  and  (d)  Hybrid  for 

10000-keV  Electrons  in  Silicon. 
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4.3.ENERGY-LOSS  SPECTRA 


In  this  section,  we  present  the  energy-loss  spectra  for  electrons  that  are  either  transmitted  or 
reflected.  The  energy-loss  spectra  discussed  below  were  obtained  for  thin  slab  problems 
according  to  the  schematic  in  fig.  1  with  exception  that  we  are  now  interested  in  energy 
dependence  rather  than  angular  dependence.  If  a  particle  is  transmitted  or  reflected  the  quantity 
tallied  is  the  initial  energy  of  the  particle  less  the  final  energy  of  the  particle  or  the  total  energy 
deposited  in  the  slab  by  the  particle. 

We  begin  with  the  energy-loss  spectra  for  100-keV  source  particles.  As  a  function  of  atomic 
number,  the  transmitted  energy-loss  spectra  do  not  change  significantly,  so  only  gold  results  are 
presented.  At  lower  energies,  it  may  not  be  feasible  to  use  a  hybrid  model  because  even  with 

2^  =  %  there  is  still  significant  disagreement  for  small  energy  transfers.  If  the  smooth  total  cross 

section  is  pushed  to  roughly  the  analog  total  cross  section,  efficiency  gains  will  be  negligible. 
Therefore,  a  mixture  of  GBFP  physics  and  analog  physics  were  tested. 
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Figure  12:  Reflected  Energy-Loss  Spectra  (a)  and  Relative  Difference  (b)  for  100-keV 

Electrons  in  Gold. 


As  seen  in  figs.  12a,  12b,  13a,  and  13a,  this  approximation  gives  good  agreement  using  only  one 
discrete  angle.  The  model  using  a  single-angle  DCS  tends  to  underestimate  the  energy-loss 
spectra.  Increasing  the  discrete  elastic  DCS  to  four  angles  results  in  very  good  agreement  with 
the  benchmark.  One  cause  of  this  could  be  the  relative  size  of  the  total  cross  section  of  elastic  to 
inelastic  scattering.  That  is,  with  a  single-angle  DCS,  elastic  scattering  is  no  longer  a  dominant 
process.  This  results  in  more  energy-loss  collisions,  so  there  are  fewer  particles  transmitted 
because  too  many  are  deposited  in  the  slab.  Hence,  the  relative  difference  plots  shows  the 
spectrum  is  underestimated  by  about  5%.  With  the  four-angle  DCS,  elastic  scattering  is  better 
represented  and  the  solution  is  improved  greatly.  In  both  cases,  there  are  very  few  particles 

contributing  beyond  Q>  so  there  are  large  oscillating  relative  differences  for  Q>\. 

Insufficient  sampling  rather  than  issues  with  the  approximation  causes  the  large  oscillations. 
These  results  were  not  included  in  the  relative  difference  plots. 
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(a)  Energy-Loss  Spectra 


A  E  (McV) 

(b)  Relative  Difference 


Figure  13:  Transmitted  Energy-Loss  Spectra  (a)  and  Relative  Difference  (b)  for  100-keV 

Electrons  in  Gold. 


For  10000-keV  source  electrons,  results  are  given  in  figs.  14a  and  14b.  We  only  present  the 
transmitted  energy-loss  spectra  because  the  reflected  energy-loss  spectra  were  overwhelmed  with 
statistical  uncertainty.  For  higher  energy  electrons,  the  elastic  DCS  has  a  greater  impact  on 
accuracy.  In  this  case,  the  single  angle  model  is  very  close  to  unity,  so  a  majority  of  the  particles 
are  pushed  through  the  slab  with  small  losses  in  energy.  Hence,  the  small  energy  losses  are 
overestimated  for  the  single-angle  approximations.  Refining  the  angular  model  improves  the 
estimation  of  the  small  energy  losses.  Additionally,  hybrid  models  are  effective  at  higher 
energies. 


(a)  Energy-Loss  Spectra  (b)  Relative  Difference 

Figure  14:  Transmitted  Energy-Loss  Spectra  (a)  and  Relative  Difference  (b)  for  10000-keV 

Electrons  in  Gold. 
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4.4.  DOSE-DEPTH  CURVES 


In  this  section,  we  present  the  dose-depth  curves  for  electrons  in  gold.  The  dose-depth  curves 
discussed  below  were  obtained  for  thin  slab  problems  according  to  the  schematic  in  fig.  1  with 
exception  that  we  are  now  interested  in  the  spatial  distribution  of  dose  deposited  throughout  the 
slab.  A  key  distinction  between  electrons  and  heavy  charged  particles  is  the  Bragg  peak  or  lack 
thereof.  Since  electrons  more  often  undergo  large-angle  scattering,  the  Bragg  peak  is  diffused 
out.  Figures  15  and  16  present  the  dose-depth  curves  for  100-keV  and  10000-keV  electrons  in 
gold. 


Depth  (cm)  x  10 4 

Figure  15:  Analog  Dose  Profile  for  100-keV  Electrons  in  Gold  Slab  with  a  Thickness  of  One 

Step. 


In  figs.  15  and  16,  there  is  no  Bragg  peak.  The  only  major  distinction  between  the  low-energy 
and  high-energy  problems  is  that  the  peak  of  the  distribution  is  pushed  further  into  the  slab 
because  high-energy  electrons  are  significantly  more  likely  to  scatter  in  forward  directions. 


33 

Approved  for  public  release;  distribution  is  unlimited. 


Depth  (cm) 

Figure  16:  Analog  Dose  Profile  for  10000-keV  Electrons  in  Gold  Slab  with  a  Thickness  of 

One  Step. 

In  figs.  17a  and  17b,  comparisons  of  the  approximations  with  the  benchmark  for  100-keV  source 
electrons  are  presented.  Results  for  four  different  approximations  are  shown  to  capture  the  effect 
of  dependence  on  angle  and  energy.  For  both  single-angle  approximations,  the  dose  is 
overestimated  in  the  first  few  cells  and  then  underestimated  throughout  the  remaining  cells. 
Regardless,  both  single-angle  approximations  are  within  roughly  3%  of  the  analog  solution. 
Further,  improvement  in  accuracy  (within  1%  of  analog)  is  achieved  by  using  a  hybrid  elastic 
scattering  DCS  or  a  mixture  of  a  four-angle  discrete  DCS  with  an  analog  inelastic  scattering 
DCS.  The  latter  being  more  efficient. 


(a)  Dose  (b)  Relative  Difference 

Figure  17:  Dose  Profiles  (a)  and  Relative  Difference  (b)  for  100-keV  Electrons  in  Gold. 
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In  figs.  18a  and  18b,  comparisons  of  the  approximations  with  the  benchmark  for  10000-keV 
source  electrons  are  presented.  Again,  for  both  single-angle  approximations,  the  dose  is 
overestimated  in  the  first  few  cells  and  then  underestimated  throughout  the  remaining  cells.  For 
the  high-energy  problem,  the  impact  of  inelastic  scattering  is  slightly  more  significant.  That  is,  it 
does  not  matter  how  accurate  the  elastic  scattering  DCS  is  if  the  inelastic  scattering  DCS  is  a 
single-energy  point.  In  fact,  the  elastic  scattering  approximation  can  be  relaxed  to  a  four-angle 
DCS  if  an  analog  inelastic  scattering  DCS  is  used.  From  an  efficiency  standpoint,  it  is 
advantageous  to  use  the  four-angle  discrete  DCS  with  an  analog  inelastic  scattering  DCS  rather 
than  a  hybrid  angle,  discrete  energy  approximation  as  it  is  1.5-20  times  more  efficient  depending 
on  the  hybrid  DCS  used. 


Figure  18:  Dose  Profiles  (a)  and  Relative  Difference  (b)  for  10000-keV  Electrons  in  Gold. 
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4.5.  TWO-DIMENSION  DOSE  DEPOSITION 


In  this  section,  two-dimension  dose  deposition  results  are  presented.  The  problem  setup  is  shown 
in  fig.  19.  For  this  simulation,  250-keV  electrons  are  normal  on  a  cube  of  silicon  surrounding  a 
small  gold  region.  The  electrons  are  normally  incident  on  the  bottom  face  of  the  cube  and 
positioned  at  the  interface  of  the  two  materials.  There  is  a  thin  layer  of  silicon  before  the 
interface  to  ensure  that  the  beam  is  spread  slightly  before  reaching  the  interface.  The  dose  is 
averaged  over  the  x-dimension. 


The  pencil  beam  is  normally  incident 
at  the  midpoint  of  the  bottom  of  the  block. 


Figure  19:  Two-Dimension  Dose  Deposition  Problem  Schematic. 

The  purpose  of  simulating  this  problem  was  to  demonstrate  the  boundary  crossing  issues 
inherent  to  condensed  history  methods.  20  presents  the  analog  benchmark.  The  actual  structure 
of  the  gold  region  imbedded  in  the  silicon  shows  up  in  the  figure.  Upon  entering  the  gold  region, 
electrons  lose  energy  rapidly.  Some  electrons  are  reflected  out  of  the  gold  and  others  never  enter 
the  gold  region.  These  electrons  that  deposit  energy  in  the  silicon  region  do  so  at  a  slower  rate 
and  the  benchmark  shows  a  smoother  gradient  in  the  dose  in  the  silicon  region. 


r 

L  * 


113  4  5  6  7 


*  to 

Dose  (keV/g) 


Figure  20:  Dose  Deposition  Benchmark  for  250-keV  Electrons  in  Silicon  Cube  Gold  with 

Insert. 
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In  fig.  21,  a  comparison  of  the  GBFP  method  and  the  Geant4  default  physics  with  the  benchmark 
is  presented.  The  two  left  plots  are  the  two  most  efficient  approximations.  The  GBFP 
approximation  shows  some  clear  discrete  artifacts  that  are  not  present  in  the  one-dimension 
results.  Since  the  single-angle  GBFP  approximation  only  allows  forward  scattering,  most  of  the 
dose  is  pushed  forward  into  the  cube  and  does  not  spread  significantly  about  the  beam-line. 
However,  the  single-angle  GBFP  approximation  does  provide  better  results  at  the  interface  than 
the  Geant4  default  physics. 


(a)  1  Discrete  Angle, 

1  Discrete  Energy 
Speed-up~35 


(b)  8  Discrete  Angles, 
1  Discrete  Energy 
Speed-up~7 


(a)  Geant4  default  physics 
with  0.04  range  factor 
Speed-up~32 


(b)  Geant4  default  physics 
with  0.004  range  factor 
Speed-up~7 


Figure  21:  Dose  Deposition  for  250-keV  Electrons  in  Silicon  Cube  With  Gold  Insert  Using 
(A)  GBFP  1  Discrete  Angle,  1  Discrete  Energy,  (b)  Geant4  Default  Physics  with  0.04  range 
Factor,  (c)  GBFP  8  Discrete  Angles,  1  Discrete  Energy,  and  (d)  Geant4  Default  Physics 

with  0.004  range  factor. 
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The  improved  GBFP  approximation  is  significantly  more  accurate  and  only  disagrees  with  the 
benchmark  very  close  to  the  source.  Again,  the  comparable  Geant4  default  physics  model  is  less 
accurate  than  the  GBFP  method  especially  near  the  interface. 
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4.6.  REFLECTION  AND  TRANSMISSION  FRACTIONS 

In  this  section,  we  present  the  reflection  and  transmission  fractions.  The  reflection  and 
transmission  fractions  discussed  below  were  obtained  for  thin  slab  problems  according  to  the 
schematic  in  fig.  1  with  exception  that  we  are  now  interested  in  the  fraction  of  particles 
transmitted  or  reflected.  For  these  calculations,  electrons  with  energies  below  5-keV  are  assumed 
to  deposit  all  energy  locally.  This  assumption  leads  to  low  reflection  results  for  low  energy 
electrons.  Though  this  doesn’t  model  reality,  this  assumption  was  consistent  between  the 
benchmark  and  the  approximations. 

In  figs.  22a,  22b,  22c,  and  22d,  the  percent  of  particles  reflected  for  the  benchmark  and  the 
relative  differences  for  three  approximations  are  presented.  The  single-angle,  single-energy 
approximation  is  within  roughly  2%  agreement  with  the  benchmark  except  for  10-keV  source 
particles.  For  10-keV  particles,  the  single-angle  is  nearly  90°,  so  after  only  two  collisions  the 
particle  could  be  reflected.  As  the  source  energy  gets  higher,  50-keV,  the  single  angles  are 
already  significantly  more  forward  peaked,  so  they  must  suffer  several  collisions  before  being 
reflected.  With  the  mixed  GBFP  DCS  and  analog  DCS  approximations,  the  10-keV  simulation 
gets  worse.  Since  the  inelastic  scattering  is  the  dominant  process  as  the  analog  inelastic  MFP  is 
simulated,  particles  likely  undergo  several  inelastic  collisions  before  experiencing  an  angular 
deflection.  As  the  particles  energy  drops  below  10-keV,  the  single-angle  DCS  is  represented  by 
an  angle  greater  than  90°  or  a  backscatter.  After  backscattering,  the  particle  undergoes  several 
more  inelastic  collisions  and  then  escapes  the  slab.  The  fact  that  the  approximation  gets  worse 
going  from  a  single-energy  inelastic  DCS  to  an  analog  inelastic  DCS,  indicates  that  there  is  some 
cancellation  of  error  in  the  fully  discrete  approximation  at  10-keV.  The  most  accurate 
approximation  uses  four  angles  and  has  very  good  agreement  with  the  benchmark. 
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Figure  22:  The  Fraction  of  Particles  Reflected  (a)  and  the  relative  difference  for  (b)  Single- 
Angle,  Single-Energy  DCS,  (c)  Single-Angle,  Analog  Inelastic  DCS,  and  (d)  Four-Angle, 

Analog  Inelastic  DCS. 
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In  figs.  23a,  23b,  23c,  and  23d  the  percent  of  particles  transmitted  for  the  benchmark  and  the 
relative  differences  for  three  approximations  are  presented.  The  single-angle,  single-energy 
approximation  is  within  roughly  1-4%  agreement  with  the  benchmark.  The  mixed  GBFP  DCS 
and  analog  DCS  approximations  are  an  improvement  over  the  fully  discrete  approximation. 
Again,  the  most  accurate  four-angle  approximation  has  very  good  agreement  with  the 
benchmark. 
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Figure  23:  The  Fraction  of  Particles  Transmitted  (a)  and  the  relative  difference  for  (b) 
Single-Angle,  Single-Energy  DCS,  (c)  Single-Angle,  Analog  Inelastic  DCS,  and  (d)  Four- 

Angle,  Analog  Inelastic  DCS. 
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5.  CONCLUSIONS 


In  many  of  the  test  cases,  the  GBFP  method  demonstrated  superb  efficiency  gains  over  analog 
with  acceptable  levels  of  accuracy.  In  comparison  with  the  Geant4  default  physics,  the  GBFP 
method  tended  to  give  more  accurate  results  for  the  same  efficiency  cost.  Extensive  testing  was 
completed  and  valuable  information  for  improvement  of  the  current  implementation  was  gained 
from  the  test  suite.  Corrections  were  made  to  the  test  code  to  account  for  the  programming  errors 
encountered  during  the  extensive  testing.  Though  the  findings  from  initial  test  suite  were 
promising,  several  test  cases  remain.  Additional  multilayer  problems  should  be  completed  to 
better  understand  the  accuracy  of  the  GBFP  method  at  boundaries.  In  addition,  only  normal 
incident  sources  were  thoroughly  tested,  so  several  off-normal  beam  sources  should  be  studied  to 
see  how  angle  of  incidence  impacts  particle  reflection.  Also,  the  partial-wave  DCSs  were  not 
thoroughly  studied.  Thorough  testing  of  the  DCS  data  will  provide  an  interesting  insight  into  the 
physics  and  more  practical  details  like  data  storage  issues. 

Implementation  of  the  advance  Monte  Carlo  method  for  charged  particles,  the  GBFP  method,  is 
nearly  complete.  Once  secondary  electron  production  is  incorporated,  implementation  of  the 
physics  will  be  complete.  Additional  testing  will  begin  to  determine  the  accuracy  of  the  GBFP 
method  when  calculating  quantities  more  relevant  to  secondaries  like  charge  deposition. 
Currently,  the  ability  of  GBFP  to  handle  secondaries  is  unknown,  but  presumed  to  be 
problematic.  Since  the  GBFP  inherently  does  not  preserve  the  total  cross  section,  and,  therefore, 
particle  number,  it  is  possible  that  only  a  mixture  of  GBFP  elastic  scattering  and  analog  inelastic 
scattering  will  result  in  accurate  and  efficient  charged  deposition  results. 

There  is  still  an  additional  implementation  detail  that  must  be  addressed.  That  is,  adaptive  GBFP 
physics.  Adaptive  GBFP  physics  will  provide  users  with  an  on-the-fly  optimization  of  efficiency 
and  accuracy  for  various  applications.  The  purpose  of  adaptivity  is  to  relieve  users  from  the  need 
to  determine  an  optimal  GBFP  representation  of  the  physics  for  their  problem  by  trial  and  error. 
We  saw  that  for  dose  calculations,  the  first  few  cells  closest  to  the  source  require  a  more  accurate 
representation  of  the  DCSs.  However,  as  the  beam  spreads  further  into  the  slab,  the  DCS  can  be 
relaxed  into  a  more  efficient  representation  and  still  maintain  sufficient  levels  of  accuracy. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


AFRL 

Air  Force  Research  Laboratory 

DCS 

Differential  Cross  Section 

GBFP 

Generalized  Boltzmann  Fokker-Plank 

MFP 

Mean  Free  Path 

UNM 

University  of  New  Mexico 
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